
%calibration
%estimating alp
cali=@(x)sum(grpsize.*( (ex_ante_allocation(x,updatek,ordering_list)-group).^2 ),'all'); 
x0=ones(size(group(1:end)))+0.5;
A = [];
b = [];
Aeq = [];
beq = [];
nonlcon=[];
lb= zeros(size(x0));
ub=[];
options = optimoptions(@fmincon,'Display','iter','Algorithm','sqp','MaxFunctionEvaluations',1300,'MaxIterations',100,'OptimalityTolerance',1e-5,'StepTolerance',1e-5,'UseParallel',true);
%
%
tic
[alp,distance_beta]= fmincon(cali,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
toc
% % 
alp_beta=alp;%[alp;5];
grpalp=alp;%[alp;5];
% 
% calculating r^*
worstofffun=@(type)interim_total(alp_beta,type,ordering_list,grpsize)-1;
worstoff=fzero(worstofffun,[0.001,1]);
r_ast_beta=zeros(1,ngroup);
parfor agent=1:ngroup
     comb=ordering_list{agent}{1};
     lower=ordering_list{agent}{2};
     q_agent=ordering_list{agent}{3};
     numberpossibilities=ordering_list{agent}{4};
     r_ast_beta(agent)=interim_allocation(q_agent,grpalp,worstoff,comb,lower,numberpossibilities); 
end
r_ast_vec_beta=r_ast_beta(ikc); %1*570 
r_fac_beta=r_ast_vec_beta(ia2);
alp_fac_beta=alp_beta(ikc);
alp_fac_beta=alp_fac_beta(ia2);
r_ast_beta
% % % 
% % % % % obtaining the fitted expected allocation 
% fitted_ea_pi=ex_ante_allocation(alp_pi,updatek,ordering_list);
% ea_fac_vec_pi=fitted_ea_pi(ikc);
% ea_fac_pi=ea_fac_vec_pi(ia2);
% % % %
% % % % revenue at r^*
% r_ast_rev=revenue(alp_pi,updatek,ordering_list,worstoff*ones(ngroup,1),r_ast_pi,alp_fac_pi,grpsize)
% % % % %revenue at r^sw
% sw_wt=worstoff_types(updatek,alp_pi,ordering_list,fitted_ea_pi,alp_fac_pi)
% r_sw_rev=revenue(alp_pi,updatek,ordering_list,sw_wt,fitted_ea_pi,alp_fac_pi,grpsize)
% propk_wt=worstoff_types(updatek,alp_pi,ordering_list,groupk/sum(fac_updatek),alp_fac_pi)
% r_propk_rev=revenue(alp_pi,updatek,ordering_list,propk_wt,groupk/sum(fac_updatek),alp_fac_pi,grpsize)
% % % % % 
% % % % % %%empirical revenue 
% emp_wt=zeros(size(endowment));
% emp_rev=zeros(nperiod,1);
% parfor tt=1:nperiod
%     empirical_share=min(max(endowment(tt,:),mini_demands'),fac_updatek'+mini_demands')-mini_demands';
%     emp_wt(tt,:)=worstoff_types(updatek,alp_pi,ordering_list,empirical_share,alp_fac_pi)';
%     emp_rev(tt)=revenue(alp_pi,updatek,ordering_list,emp_wt(tt,:)',empirical_share,alp_fac_pi,grpsize);
% end
% emp_rev

% test_share=zeros(1,57);
% test_share(find(fac_updatek==groupk(1),1,'first'))=1.7431-1.5691;
% %test_share(find(fac_updatek==groupk(4),2,'first'))=groupk(3);
% %test_share(find(fac_updatek==groupk(4),1,'last'))=1-grpsize(2)*groupk(2);
% test_wt=worstoff_types(updatek,alp_pi,ordering_list,test_share,alp_fac_pi)'
% test_rev=revenue(alp_pi,updatek,ordering_list,test_wt',test_share,alp_fac_pi,grpsize)
% 

% groupnumber=2;
% groupk(groupnumber)
% comb=ordering_list{groupnumber}{1};
% lower=ordering_list{groupnumber}{2};
% q_agent=ordering_list{groupnumber}{3};
% numberpossibilities=ordering_list{groupnumber}{4};
% interim_allocation(q_agent,alp_pi,0,comb,lower,numberpossibilities)
      
%%boundary solution: 
%inforents=info_rent(alp_pi,updatek,ordering_list)
%corner_share=zeros(1,57);
%corner_share(find(fac_updatek==groupk(1),1,'first'))=1;
%test_share(find(fac_updatek==groupk(4),2,'first'))=groupk(3);
%test_share(find(fac_updatek==groupk(4),1,'last'))=1-grpsize(2)*groupk(2);
%corner=worstoff_types(updatek,alp_pi,ordering_list,test_share,alp_fac_pi)';
%corner_rev=revenue(alp_pi,updatek,ordering_list,corner_wt',test_share,alp_fac_pi,grpsize)
% % calculating hybrid r^*
% effective_endow=norm2-mean(sum(endowment_adj,2));
% worstofffun=@(type)interim_total(alp_pi,type,ordering_list,grpsize)-(1-effective_endow);
% worstoff_h=fzero(worstofffun,[0.001,1]);
% r_h_pi=zeros(1,ngroup);
% parfor agent=1:ngroup
%      comb=ordering_list{agent}{1};
%      lower=ordering_list{agent}{2};
%      q_agent=ordering_list{agent}{3};
%      numberpossibilities=ordering_list{agent}{4};
%      r_h_pi(agent)=interim_allocation(q_agent,grpalp,worstoff_h,comb,lower,numberpossibilities); 
% end
% r_h_vec_pi=r_h_pi(ikc); %1*570 
% r_h_pi_fac=r_h_vec_pi(ia2);
% r_h_pi
% 
% r_h_rev=revenue(alp_pi,updatek,ordering_list,worstoff_h*ones(ngroup,1),r_h_pi,alp_fac_pi,grpsize)


% calculating hybrid r^*, exogeneously given over endowment
% overendowed_fac=find(average_endowment_fac>=fac_updatek); %agents with overendowments
% normal_fac=find(average_endowment_fac<fac_updatek); 
% for gg=1:ngroup
%     grpsize_over(gg)=grpsize(gg)-sum(ismember(fac_updatek(overendowed_fac),groupk(gg)));
% end
% effective_endow=1-(sum(average_endowment_fac(overendowed_fac))-sum(mini_demands(overendowed_fac)));
% worstoff_over=@(type)interim_total(alp_pi,type,ordering_list,grpsize_over)-effective_endow;
% worstoff_h_over=fzero(worstoff_over,[0.001,1]);
% r_h_over_pi=zeros(1,ngroup);
% parfor agent=1:ngroup
%      comb=ordering_list{agent}{1};
%      lower=ordering_list{agent}{2};
%      q_agent=ordering_list{agent}{3};
%      numberpossibilities=ordering_list{agent}{4};
%      r_h_over_pi(agent)=interim_allocation(q_agent,grpalp,worstoff_h_over,comb,lower,numberpossibilities); 
% end
% r_h_over_vec_pi=r_h_over_pi(ikc); %1*570 
% r_h_over_pi_fac=r_h_over_vec_pi(ia2);
% %r_h_over_pi
% r_h_over_pi_fac_ed=r_h_over_pi_fac;
% r_h_over_pi_fac_ed(overendowed_fac)=average_endowment_fac(overendowed_fac)-mini_demands(overendowed_fac);
% % sum(r_h_over_pi_fac_ed)
% over_wt=worstoff_types(updatek,alp_pi,ordering_list,r_h_over_pi_fac_ed,alp_fac_pi)
% r_h_over_rev=revenue(alp_pi,updatek,ordering_list,over_wt,r_h_over_pi_fac_ed,alp_fac_pi,grpsize)
% 
% 
% % 
% %trading volume, only account for the non-artifical facilities
% cap_mean_eff_trading_vol=zeros(nperiod,1);
% cap_eff_trading_vol_5=zeros(nperiod,1);
% cap_eff_trading_vol_95=zeros(nperiod,1);
% mean_eff_trading_vol=zeros(nperiod,1);
% eff_trading_vol_5=zeros(nperiod,1);
% eff_trading_vol_95=zeros(nperiod,1);
% mean_ast=zeros(nperiod,1); ast_5=zeros(nperiod,1);ast_95=zeros(nperiod,1);
% mean_sw=zeros(nperiod,1); sw_5=zeros(nperiod,1);sw_95=zeros(nperiod,1);
% for t=1:nperiod
%      alloc=endowment(t,:);%*supply_yearly.supply(t);
%      adjalloc=min(max(endowment(t,:),mini_demands'),fac_updatek'+mini_demands');
%      ast_alloc=r_fac_pi+mini_demands'; %1*57
%      sw_alloc=ea_fac_pi'+mini_demands';
%      ndraw=1000;
%      nagent=nfac;
%      vol=zeros(1,ndraw);
%      voladj=zeros(1,ndraw);
%      ast=zeros(1,ndraw);
%      sw=zeros(1,ndraw);
%      for rep=1:ndraw
%         types=zeros(nagent,1);
%         for agent=1:nagent
%             a=alp_fac_pi(agent);
%             types(agent)=invdist(rand(1),a); 
%         end
%         [~,ranking]=sort(types,'descend');
%         allocation=realized_allocation(ranking,fac_updatek)'+mini_demands';
%         vol(rep)=sum(abs(allocation-alloc))/2;%-abs(allocation(57)-alloc(57));
%         voladj(rep)=sum(abs(allocation-adjalloc))/2;%-abs(allocation(57)-adjalloc(57));
%         ast(rep)=sum(abs(allocation-ast_alloc))/2;%-abs(allocation(57)-ast_alloc(57));
%         sw(rep)=sum(abs(allocation-sw_alloc))/2;%-abs(allocation(57)-sw_alloc(57));
%      end
%      mean_eff_trading_vol(t)=mean(vol);
%      eff_trading_vol_5(t)=prctile(vol,5);
%      eff_trading_vol_95(t)=prctile(vol,95);
%      mean_ast(t)=mean(ast);
%      ast_5(t)=prctile(ast,5); 
%      ast_95(t)=prctile(ast,95);
%      mean_sw(t)=mean(sw);
%      sw_5(t)=prctile(sw,5); 
%      sw_95(t)=prctile(sw,95);
%      cap_mean_eff_trading_vol(t)=mean(voladj);
%      cap_eff_trading_vol_5(t)=prctile(voladj,5);
%      cap_eff_trading_vol_95(t)=prctile(voladj,95);
% end
% 

% % functions: 
% 

%v1
% function[meanvol,vol5,vol95,mean_ast,ast_5,ast_95]=trading_volume(k_agent,alloc,r,alp_agent) 
%     ndraw=1000;
%     nagent=length(k_agent);
%     vol=zeros(1,ndraw);
%     ast=zeros(1,ndraw);
%     alp=alp_agent;
%     for rep=1:ndraw
%         types=zeros(nagent,1);
%         for agent=1:nagent
%             a=alp(agent);
%             types(agent)=invdist(rand(1),a); 
%         end
%         [~,ranking]=sort(types,'descend');
%         allocation=realized_allocation(ranking,k_agent);
%         vol(rep)=sum(abs(allocation-alloc))/2;
%         ast(rep)=sum(abs(allocation-r))/2;
%     end
%     meanvol=mean(vol);
%     vol5=prctile(vol,5);
%     vol95=prctile(vol,95);
%     mean_ast=mean(ast);ast_5=prctile(ast,5);ast_95=prctile(ast,95);
% end
% function[rev]=revenue(grpalp,grpk,ordering_list,worstofftypes,alloc,alpfac,grpsize)
%      nagent=length(alloc);
%      ngroup=length(grpk);
%      rev=zeros(nagent,1);
%      for agent=1:nagent
%          if ngroup==nagent
%             groupnumber=agent;
%         else
%             groupnumber=find(grpalp==alpfac(agent));
%         end
%         comb=ordering_list{groupnumber}{1};
%         lower=ordering_list{groupnumber}{2};
%         q_agent=ordering_list{groupnumber}{3};
%         numberpossibilities=ordering_list{groupnumber}{4};
%         a=grpalp(groupnumber);
%         w=worstofftypes(agent);
%         fun2=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities).*buyertypetimesdensity(a,type);
%         fun1=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities).*sellertypetimesdensity(a,type);
%         rev(agent)=integral(fun1,0,w)+integral(fun2,w,1)-max(0,min(alloc(agent),grpk(groupnumber)))*w;
% %         m=10;
% %         x0=0;
% %         x1=1;
% %         dx=(x1-x0)/m;
% %         int=0;
% %         parfor i=1:m
% %             int=int+ fun(x0+(i-0.5)*dx)*dx
% %         end
% %         eq(agent)=int;
%      end
%      if ngroup==nagent
%         rev=sum(rev.*grpsize);
%      else
%         rev=sum(rev);
%      end
% end

% function[wt]=worstoff_types(groupk,grpalp,ordering_list,alloc,alpfac) % note alloc has to be the same size as grpk
%     nagent=length(alloc);
%     ngroup=length(groupk);
%     wt=zeros(nagent,1);
%     for agent=1:nagent
%         if ngroup==nagent
%             groupnumber=agent;
%         else
%             groupnumber=find(grpalp==alpfac(agent));
%         end
%         comb=ordering_list{groupnumber}{1};
%         lower=ordering_list{groupnumber}{2};
%         q_agent=ordering_list{groupnumber}{3};
%         numberpossibilities=ordering_list{groupnumber}{4};
%         fun=@(type)min(alloc(agent),groupk(groupnumber))-interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
%         if fun(0)<0
%              wt(agent)=0;
%         elseif fun(1)>0
%             wt(agent)=1;
%         else
%             wt(agent)=fzero(fun,[0,1]);
%         end
%     end
% end

function[eq]=ex_ante_allocation(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)density(a,type).*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0,1);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end
function[eq]=info_rent(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0.01,0.99);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end

function[suma]=interim_total(grpalp,type,ordering_list,grpsize) %remember to change interim total 
    nagent=length(grpalp);
    suma=0;
    parfor j=1:nagent
        comb=ordering_list{j}{1};
        lower=ordering_list{j}{2};
        q_agent=ordering_list{j}{3};
        numberpossibilities=ordering_list{j}{4};
        suma=suma+grpsize(j)*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
    end    
end

%v1:
% function[agent]=interim_allocation(grp,grpalp,type,grpk,comb,lower) %alp needs to be a col vector   
%     q_agent=agent_allocation(grp,comb,grpk);
%     alp=grpalp;
%     higher=comb;
%     ntype=length(type);
%     agent=zeros(size(type));
%     parfor tt=1:ntype
%         prob11=prob1(alp,type(tt),higher);
%         prob22=prob2(alp,type(tt),lower);
%         agent(tt)=sum(prob11.*prob22.*q_agent,'all');
%     end
% end

%v2: use v2 of ordering_list
function[agent]=interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities) %alp needs to be a col vector   
    %q_agent=agent_allocation(grp,comb,grpk);
    alp=grpalp;
    higher=comb;
    ntype=length(type);
    agent=zeros(size(type));
    parfor tt=1:ntype
        prob11=prob1(alp,type(tt),higher);
        prob22=prob2(alp,type(tt),lower);
        agent(tt)=sum(prob11'.*prob22'.*q_agent.*numberpossibilities,'all');
    end
end

function[q]=realized_allocation(order,fac_updatek)
    nagent=length(order);
    q=zeros(nagent,1);
    rest=1;
    allocated=0;
    k=fac_updatek;
    for rk=1:nagent
        agent=order(rk);
        q(agent)=min(rest,k(agent));
        allocated=allocated+q(agent);
        rest=max(0,1-allocated);
    end
end


function[z1]=density(a,type)
    z1=pdf('Beta',type,a+1,a+1);
end
% function[z1]=buyertypetimesdensity(a,type)
%     z1= a.*(type).*(1-type).^(a-1)-(1-type).^a;
% end
% 
% function[z1]=sellertypetimesdensity(a,type)
%     z1= a.*(type).*(1-type).^(a-1)+(1-(1-type).^a);
% end

% function[t]=invdist(u,a) %given quantile, output type
%     t=  1-(1-u)^(1/a);
% end
function[z2]=prob1(a,type,higher)
    z2=prod((1-cdf('Beta',type,a+1,a+1)).^higher,1);
end

function[z3]=prob2(a,type,lower)
    z3=prod((cdf('Beta',type,a+1,a+1)).^lower,1);
end


